QGP flow fluctuations and the caracteristics of higher moments 
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The dynamical development of expanding Quark-gluon Plasma (QGP) flow is studied in a 3+1D 
fluid dynamical model with a globally symmetric, initial condition. We minimize fluctuations aris- 
ing from complex dynamical processes at finite impact parameters and from fluctuating random 
initial conditions to have a conservative fluid dynamical background estimate for the statistical dis- 
tributions of the thermodynamical parameters. We also avoid a phase transition in the equation 
of state, and we let the matter supercool during the expansion. Then central Pb+Pb collisions 
at y/SNN = 2.76 TeV are studied in an almost perfect fluid dynamical model, with azimuthally 
symmetric initial state generated in a dynamical flux-tube model. The general development of ther- 
modynamical extensives are also shown for lower energies. We observe considerable deviations from 
a thermal equilibrium source, changing skewness and kurtosis by time depending on beam energy 
as a consequence of the fluid dynamical expansion arising from a least fluctuating initial state. 

PACS numbers: 12.38.Mh, 25.75.-q, 25.75.Nq, 51.20.+d 
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COLLECTIVE FLOW OF GLOBAL 
SYMMETRY 



In heavy ion collisions collective flow has been mea- 
sured and azimuthal asymmetry was determined from v\ 
to vs ■ At the highest energies in central collisions fluctua- 
tions dominate arising from fluctuating initial configura- 
tions, and the most dominant flow harmonic is v^. These 
are collective flow fluctuations and have no direct connec- 
tion to the fluctuations arising from a phase transition in 
the Equation of State (EoS) 

The fluctuations arising from the pure fluid dynam- 
ics without hadronization are studied in this work. We 
present a set of calculations, what kind of effects a possi- 
bly most conservative relativistic fluid dynamical model 
exhibits in higher moments of statistical parameters for 
extensive densities, as skewness and kurtosis, without in- 
cluding any effect associated with a phases transition, in 
the EoS, in the transport properties, or in special ther- 
modynamical phase space trajectories 043; or m special 
freeze out mechanisms. This study complements a large 
number of studies with the opposite goal, aiming to an- 
alyze the consequences of the above mentioned effects. 
Contrary to the effort to eliminate all these effects, we 
still obtain energy dependent changes in the skewness 
and kurtosis. 

We use an EoS without a phase transition and include 
the possibility of supercooling: we choose the MIT Bag 
model with parameters fixed to the initial values (two 
flavors and massless quarks and gluons, and the bag con- 
stant is B = 0.397 GeV/Jm 3 ). In addition to omitting 
the freeze out and its effects, we also avoid to take into 
account viscosity (except the inavoidable numerical vis- 
cosity) , thus also the temperature dependence of viscosity 
near the critical point [8|, which may lead to additional 



changes of the critical fluctuations in a viscous fluid dy- 
namical evolution. 

In peripheral collisions, on the other hand, the ini- 
tial asymmetry is dominated by the almond shape of the 
participant matter. This results in a strong elliptic flow, 
while the directed flow measured in the ALICE TPC ap- 
peared to be weak and dominated by random fluctua- 
tions. 

Computational Fluid Dynamics (CFD) predictions in- 
dicated a new directed flow structure: due to the large 
angular momentum of the initial state in peripheral colli- 
sions the anti-flow peak observed at high SPS and RHIC 
energies is rotating forward, and at sufficiently high beam 
energy, v\(y) will start to peak at positive rapidities, i.e. 
on the same side where the projectile spectator residues 
arrive after the collision This happens because the 
initial angular momentum leads to a faster rotating initial 
system, and this rotation moves the dominant directed 
flow peak forward before the expansion from the pressure 
would slow down the rotation. The observation of this 
peak is not easy because of the beam directed fluctua- 
tions of the initial state. 

At high energies the dimensionless shear viscosity over 
the entropy density, ij/s, of the QGP is becoming small 
Q5, and ri/s as a function of temperature has a mini- 
mum at the critical temperature [111 ]. So the Reynolds 
number may exceed one, and turbulent phenomena may 
start to occur. On the other hand, ij/s at a critical point 
does not necessarily have a minimum, since the dynami- 
cal universality class of a possible critical point of QCD is 
the H-model in Hohenberg and Halperin's classification 
[l~2l [l3| . the shear viscosity may diverge at a possible 
QCD critical point. This would damp instabilities. 

Recently in the same CFD model with the Particle in 
Cell (PIC) method, it was observed that in peripheral col- 
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lisions with low viscosity (and low numerical viscosity) a 
Kelvin- Helmholtz (KH) instability starts to develop [Tj] , 
which enhances the rotation effect and the spatial vari- 
ance of the flow pattern. The flow effects depend on the 
initial state profile and on the viscosity. Turbulence ap- 
pears only for small viscosity, which indicates the critical 
point of the matter fll| , and it is a sensitive measure of 
viscosity and its minimum at the critical point. 

Our CFD simulations of the LHC heavy ion collisions 
suggest that collective directed v\{y) flow function can 
be measured if the Globally asymmetric flow component 
and the random flow arising from the initial state ycM- 
fluctuations can be separated. In hydrodynamical calcu- 
lations we see that the v\ Global flow can change the peak 
position to "forward" with increasing beam energy and 
initial angular momentum. This is a result of our tilted 
initial state with shear flow fi~5l . , in which the angular 
momentum from the increasing beam momentum may 
supersede the expansion driven by the pressure. 

The above described phenomena contribute to an in- 
creased spatial spread of the matter during the collision. 
These effects are present even if we do not have a phase 
transition in our EoS [17|. Our present goal is to de- 
termine the lowest possible deviation from an ideal ther- 
mal source in a least fluctuating CFD evolution. Thus, 
we eliminate random initial state fluctuations and all az- 
imuthal asymmetries in a head-on collisions to obtain a 
most symmetric distribution, and study the spatial fluc- 
tuation of thcrmodynamical quantities in such a system. 
This then can be compared to the effects caused by the 
phase transition QH^. 

Most random fluctuations lead to close to Gaussian 
distribution, nevertheless the dominant fluid dynamical 
expansion, even if all special sources are eliminated or 
minimized may lead to more complex non-gaussian fluc- 
tuations and higher statistical moments. 

For a realistic reaction model we have to describe the 
final stage of the reaction also. We have a Multi Module 
Model approach to describe high energy heavy ion colli- 
sions in the RHIC and LHC energy range. Then from the 
locally equilibrated QGP we have to form hadrons. We 
do not assume that the hadronization happens in chem- 
ical equilibrium as this would take too long time [2l|, [22[ 
and would not allow for baryons of high strangeness. 
Thus we use the simplest bag model approach, which in 
the pure QGP domain, yields similar results to more de- 
tailed parametrizations fitted to lattice predictions [23| . 

In order to be able to hadronize rapidly we have to 
assume a fast, non-equilibrium hadronization and freeze 
out with either a Cooper-Frye based method with a lo- 
cal sudden change or a sudden transition to a parton 
and hadron cascade model, e.g. 24], which can describe 
rapid hadronization without the assumption of chemi- 
cal equilibration. This final stage should additionally 
increase the deviations from the local statistical equilib- 
rium. Thus, the final observed (or calculated) particle 
distribution, which contains already the influence of the 
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FIG. 1. (color online) The trajectory of fluid dynamical de- 
velopment of QGP fluid at different beam energies as indi- 
cated in the figure. Open squares belong to central Pb+Pb 
collisions, the others are central Au+Au collisions, these are 
calculated at a cell size resolution of dx — dy = dz — 0.575 
frn, and time step At — 0.04233 fm/c. The hadronic freeze 
out curve [25[ is indicated by a full black line. The CFD evo- 
lution is calculated well beyond this curve. This is possible 
as the CFD model can describe supercooled QGP fluid also. 
The viscosity is minimal and only the numerical viscosity is 
considered in the calculations. These can be performed down 
to FAIR and NICA energies, although the use of supercooled 
QGP EoS has constrained validity at these low energies. 



rapid hadronization and freeze out, should be compared 
to the basic fluctuations arising from a (least fluctuating) 
CFD distribution estimate. 

It is also important to mention that different thermo- 
dynamical parameters (especially intensives and exten- 
sives) do not have to show the same critical fluctuation 
properties, so we have to study the fluctuations of sev- 
eral parameters. Furthermore the statistical physics esti- 
mates assume a single thermal source at or near the criti- 
cal point, while we estimate here also the effects of spatial 
fluctuations, which arise from a dynamically expanding 
fluid flow even in the least fluctuating configuration. 



II. FLOW DEVELOPMENT 

The dynamically developing flow pattern leads to a 
spatial distribution of all thermodynamical quantities, 
while the system expands rapidly. We assume that the 
most probable scenario is a pre-equilibrium development 
described by linear flux-tube expansion independently at 
each point of the transverse, [a;, y], plane until local equi- 
librium is reached at a space-time hypersurface. By this 
time at high energies we reach a (nearly) equilibrated 
Quark-gluon Plasma (QGP) state, which then expands 
and supercools. This intermediate stage is described with 
a CFD model using the PIC method. Finally the super- 
cooled QGP can hadronize rapidly and almost simulta- 
neously it freezes out. This final stage of the reaction can 
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FIG. 2. (color online) The average specific energy density 
in the calculational frame (e CF ), the average specific internal 
energy density in the local rest frame (s LR ) and the average 



specific kinetic energy (e = e 



CF 



e ) are calculated by 



the PIC hydro model in central collisions at ^/snn = 2.76 
TeV. The star indicates eo = 0.938 GeV, the initial specific 
internal energy before collision. The error bars indicate the 
variance, a. 
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FIG. 3. (color online) The average specific energy density 
in the calculation frame (e CF ), the average specific internal 
energy density in the local rest frame (s ) and the average 
specific kinetic energy (e K ) are calculated by the PIC hydro 
model in central collisions at -^/sjvjv = 200 GeV. eo = 0.938 
GeV is the initial internal energy before collision. The kinetic 
energy is the difference between e CF and e LR . 



be described by a non-equilibrium model. 

The CFD model is using N fluid cells, i = 1,2, ...,N, 
where with time and expansion the number of fluid cells 
is increasing. Thus, each cell carries less and less baryon 
charge with time, and a different amount. Therefore to 
calculate the volume average and the distribution over 
the volume we weight the fluid cells by the amount of 
baryon charge they carry. 

Thus the weighted average of a quantity x is defined 
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FIG. 4. (color online) The average specific energy density 
in the calculation frame (e CF ), the average specific internal 
energy density in the local rest frame (e LR ) and the average 
specific kinetic energy (e K ) are calculated by the PIC hydro 
model in central collisions at ^/sjvjv = 39 GeV. eo = 0.938 
GeV is the initial internal energy before collision. The kinetic 
energy is the difference between e CF and e LR . Observe that at 
7 fm/c the specific internal energy and kinetic energy change 
stalls indicating that the majority of the cells reached zero 
pressure, the FO boundary. This coincides with the crossing 
point of the FO curve in Fig. [TJ indicating that beyond that 
point the assumed development in the supercooled QGP is 
overstretching the applicability of the CFD calculation. 



as 
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where 
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and Vi is the volume of i th fluid cell, nf F is the baryon 
density in the Calculational Frame (CF), which is the 
cm. frame in the present calculations, and nf F = nf R ji 
(where nf R is the baryon density in Local Rest Frame), 
so that N to t = J2i n i R Ji v n therefore J2i w» = 1. 

The CFD stage of this development is shown in Fig. [TJ 
where the time development of the average thermody- 
namical quantities is shown in the temperature, T and 
baryon chemical potential, n, plane and the full black 
line is the hadronic freeze out curve. In the present work 
we analyze the CFD stage, and analyze statistically the 
space-time development of both the intensive and specific 
extensive thermodynamical variables. We study central 
collisions only, to avoid the effects from azimuthal flow 
asymmetries and from particle emission from projectile 
and target residues (spectator evaporation) [26j |. 

For a variable x the averages and various order mo- 
ments of cell-by-cell distributions can be written as 
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where P(x) is the spatial distribution weighted by the 
baryon charge density in the CF. The spatial variance, 
the skewness and the kurtosis can be obtained from these 
moments: 
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By using this average, first we can calculate specific ex- 
tensives, which are governed by strict conservation laws. 
The total baryon charge, energy and momentum conser- 
vations are governed by the continuity equation and by 
the relativistic Euler equation. 
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FIG. 5. (color online) The average baryon charge density 
in the calculation frame (n CF ) and in the local rest frame 
{n LR ) calculated by the PIC hydro model in central Pb+Pb 



2.76 TeV. The error bars indicate the 
variance, er, for n oi< . The variance is about the same for the 
invariant specific density, n LR . 



collisions at ^/snn - 
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and as a consequence, the total momentum in the cen- 
ter of mass (cm.) frame should remain zero during the 
development, while the average specific energy per net 
nucleoli number 



too 

< e >= —— q- = const. 



(9) 



should remain constant in CF. 

In most experiments the observable is the total charged 
particle multiplicity. This number is proportional to the 
energy density and not the net baryon charge. That is 
the reason we study the energy density. The energy is 
characterized by the "specific energy density", i.e., by 
the energy per unit net baryon charge. As you can see, 
the different parts of the energy development undergo 
significant changes with time, which may contribute to 
the statistical properties of the produced particles. 

The average specific energy density in CF can be ex- 
pressed as: 
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Similarly we can obtain the average specific energy den- 
sity in the local rest frame (LR): 



< e LR >- 
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These quantities are shown in Figs. I2I3I4I Most im- 
portantly the average specific energy remains constant 
during the time development, and its value is the initial 
beam energy, < e CF >— ^^/sjvat = 1-38 TeV. Due to 
the CFD expansion, the average specific internal energy, 



£ lr _ e j n ^ decreases as the system expands, while the 
expansion leads to increased average specific kinetic en- 
ergy, e K = e CF — e LR . Of course these quantities vary 
in the space-time during the CFD evolution. Their spa- 
tial variances, a x , are shown in eq. (U), is also indicated 
by the error bars in Figs. I2I3I4I Although the average 
of e CF remains constant its variance is increasing due to 
the expansion, which generates increasing number of low 
density cells. The lower energies give similar results. The 
three figures demonstrate how much part of the energy 
is converted into kinetic energy. At the highest energy 
by the time of 8 fm/c, half of the total available is con- 
verted into flow while at lowest shown beam energy, at 
yfsWjj = 39 GeV, it is only 20%. 

Similarly, the net baryon number, N tot = Y^i n \ n iNii 
is exactly conserved in the calculations, as the marker 
particles carry fixed baryon charge and these are con- 
served until they are in the calculation grid. At the same 
time the average baryon charge density is decreasing. We 
can characterize this by the invariant scalar, LR, baryon 
density as well as by the CF baryon density. Their aver- 
ages can be expressed as: 



<n >=2^ n l 



-CF 



Vi 



N t , 



EC" 



CF\2 



Nto; 



(12) 



Similarly we can obtain the average of the invariant scalar 
baryon density 



< n LR >-- 



N t , 
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liVi 

Ntot 
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The time dependencies of the two densities are shown 
in Fig. [5j We observe that both densities decrease with 
time, and the density calculated in CF is larger than the 
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FIG. 6. (color online) The time development of the skew- 
ness and kurtosis of the distribution of the net baryon density, 
n CF , n LR , for central Pb+Pb collision at yi^T = 2.76 TeV. 
The skewness of the baryon charge is positive for both densi- 
ties indicating an increasingly longer high density tail of the 
distribution. The kurtosis is initially negative, and turns pos- 
itive only at very late stages. This is arising from the initial 
condition where the baryon charge is uniformly distributed 
in each longitudinal "streak". Around 4 fm/c n CF is visibly 
smaller than n LR indicating that the apparent density in the 
CF is more uniform than the invariant scalar density. 
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FIG. 7. (color online) The time development of the kur- 
tosis, K, and skewness, S, of the distribution of the specific 
energy density in the CF, e CF , for central Pb+Pb collisions 
at y/SNN = 2.76 TeV. The skewness is positive, decreasing 
towards zero with time, indicating a longer high energy tail, 
especially at early stages where we have more cells with high 
energy QGP, and high kinetic energy as discussed in ref. [19l ]. 
The kurtosis is initially positive indicating a high energy den- 
sity QGP with small spread but decreases with time and it is 
becoming negative at 2 fm/c. This is caused by populating 
low energy density supercooled states at later stages. 



density calculated in the LR. Their variance is decreasing 
with time. 



A. Late stages of expansion 

As mentioned above the EoS of our perfect fluid dy- 
namical model is ideal QGP, in the form of the MIT Bag 
model. The parameter of the Bag constant is fixed until 
the local pressure is positive. In order to be able to cal- 
culate the continued expansion in the supercooled QGP 
we set and fix the pressure to zero and decrease the Bag 
constant so that the expansion in the supercooled state 
remains adiabatic 



27] 



Usually this happens only in very few cells before the 
estimated average Freeze Out (FO) time (less than 10% 
of the cells), but we did continue the CFD calculations 
well into the supercooled state when the zero-pressure 
cells amounted to 30-40% of the total volume. 

In some of the thermodynamical variables this change 
is exhibited by a change of the development trend line. 
At the same time the energetic characteristics of the EoS 
are realistic, so the overall development and the basic 
quantities are well estimated even during the supercooled 
stage. 



III. SKEWNESS AND KURTOSIS 

In the present work we assume a single QGP phase. 
We study the features the energy density and net baryon 
density variations exclusively arising from the fluid dy- 
namics. This dynamical change is not observable di- 
rectly, only at the freeze out hypersurface which is si- 
multaneous or close to the hadronization. At this point 
the hadronization may significantly modify the statisti- 
cal fluctuations [20] . This work evaluates the fluctuations 
from the previous fluid dynamical evolution. In this part 
we study the skewness and kurtosis of the specific energy 
density and the baryon density according to Eq. ([5]) and 
Eq. ([6]). In Fig. [6l we can see that the kurtosis of the 
baryon density calculated in different frames are simi- 
lar, the kurtosis is negative at first and then turns to be 
positive around t = 13 fm/c. The skewness is allover pos- 
itive. At t = 4 fm/c, both densities have a minimum, but 
the skewness value in the local rest frame is almost twice 
larger than in the calculational frame. This indicates 
that the contribution of the flow makes the distribution 
close to Gaussian. 

In Fig. [7] we can see that the kurtosis of the specific 
energy in CF, e , changes sign from positive to nega- 
tive, while the skewness is always greater than zero. This 
is not the same when compared to the results calculated 
in the local rest frame, which are shown in Fig. [8j where 
skewness changes sign from negative to positive. 

The change of statistics with the changes of the EoS is 



6 



4 



2 



-4 



1 1 1 1 1 


1 1 1 


■ i 1 i 1 i 1 


£ LR 
o S 




A A A 
A A 

A 

AA A " 


A K 




A 

^^^^ ' 







A A^ 




^A^4 





~1 1 1 ' 1 ' 1 ' 1 ' 1 1 1 1 1 1 

2 4 6 8 10 12 14 16 



t (fm/c) 

FIG. 8. (color online) The time development of the kurto- 
sis, K, and skewness, S, as in Fig. [7] but for the invariant 
scalar specific energy density, e LR , distribution. In case of 
the invariant scalar specific energy distribution, the skewness 
changes sign at 7 fm/c. Unlike in Fig. [7] the skewness is 
initially negative indicating a wider spread in the invariant 
energy density distribution, which can also be attributed to 
the presence of QGP. 



best seen in Fig. [5] for the invariant specific energy dis- 
tribution e LR , although the amplitude of the change is 
small. The quantity e LR , does not include the contribu- 
tion from the kinetic energy, thus it is the best measure 
of changes in the EoS. The changes are still observable 
contrary to the fact that the QGP to HM transition is 
not included in the EoS, but we still have a transition 
from ideal QGP to the supercooled, zero pressure QGP, 
where the energy density decreases with the decreasing 
bag constant. 

It is interesting to compare the skewness of the two 
types of energy densities in Fig. [9l e CF and e LR . While 
e CF , which contains all kinetic energy, has always a pos- 
itive skewness, i.e. the distribution is more uniform, the 
internal energy, e LR , which increases gradually, has ini- 
tially negative skewness, meaning a more spread out, 
fluctuating distribution, which becomes positive later, ac- 
tually at and after the estimated FO time, when the mat- 
ter is supercooled and characterized by zero pressure and 
reduced Bag constant. We can also note that the cross 
point of the skewness in CF and in LR is around 8 fm/c, 
which is the same as in Fig. f2l where the trajectories 
of the Local Rest energy and the Kinetic energy cross. 
The experimental observations do not show clearly the 
expected phase transition behavior [20| . That is why we 
study alternative origins of the fluctuations. The skew- 
ness and kurtosis changes are of the same order of mag- 
nitude as from critical fluctuations, therefore both effects 
should be studied. 

In a multi-module or hybrid-model construction (e.g 
where the PIC hydro stage is matched to a parton and 
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FIG. 9. (color online) The time development of the skewness, 
K, for the CF specific energy density, e LR , and the invariant 
scalar specific energy density, e LR , distribution. The invariant 
scalar specific energy distribution changes sign at the typical 
FO moment, 8 fm/c. 



hadron cascade model PACIAE), the flow features are 
matched (28j to a subsequent dynamical model which de- 
scribes dynamical, non-equilibrium, rapid hadronization. 
These types of models can describe realistically the statis- 
tical properties and a dynamical phase transition, which 
provide the hadron distribution in the final stage. This 
stage would then explicitly describe the random fluctua- 
tions arising from the phase transition also. 

In ref. [29j a mixed particle method is introduced, 
which could separate the fluctuations arising from local 
critical fluctuations. The mixed events are actually elim- 
inating two particle correlations, and only the single par- 
ticle distributions remain. Thus for central events these 
are mainly local correlations which may arise from local 
fluctuations caused by energy and baryon charge clus- 
tering in a phase transition. The method separates the 
consequences of such correlations. 

This method can be used both in hybrid model calcu- 
lations and in experiments, to separate the fluctuation 
effects from the collective flow and the phase transition 
dynamics. 

IV. CONCLUSIONS 

We studied in a fluid dynamical model the time de- 
velopment of spatial distribution averages and variances 
of thermodynamical variables in high energy heavy ion 
collisions. We assumed an Equation of State for Quark 
Gluon Plasma only, but including supercooled QGP also 
where the pressure is dropped to zero. We studied cen- 
tral collisions, without initial state fluctuations to min- 
imize fluctuations arising from complex anisotropic flow 
patterns. Including the supercooled QGP we could de- 
scribe the late stages of the collisions expanding up to 
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the phase transition boundary [25[ and beyond (Fig. [TJ. 
In addition to the average temperature and chemical po- 
tential we also studied the developments of average en- 
ergy densities and baryon densities as well as their vari- 
ances (Figs. 1213 1415 P - The dynamical developments of 
these variables showed the expected, monotonic dynam- 
ical behavior, even beyond the physical FO times, where 
we overstretched the applicability of the fluid dynamical 
model. 

Our FD model did not include the rapid hadroniza- 
tion and the random generation of hadrons, which would 
generate critical fluctuations in the vicinity of the critical 
point of the phase transition. 

Interestingly the higher statistical moments, the Kur- 
tosis and Skcwncss still showed a non monotonic behavior 
(Figs. [51 El [HIES)- We presented these for the highest en- 
ergy collisions, where the applicability of the applied EoS 
is the least questionable. These higher moments of spe- 
cific energy densities show changing signs of the Kurtosis 
and Skewness. 

Apart of the fact that the FD model provides a spatial 
distribution, Fig. [5] shows that the variation of higher 
moments is also sensitive to the FO time. This time can- 



not be securely determined from within the FD model, 
and the final hadronization and freeze-out should be de- 
scribed by the last model stage of the hybrid model. 
Thus, we did not study the excitation function of the 
higher moments, because an arbitrary choice of the FO 
time may result in different results. On the other hand 
the time dependence of higher moments obtained here 
(Figs. 6-9) is different from the straightforward expecta- 
tion arising from a dynamical phase transition [20(. 

The spatial fluctuations of specific energy and baryon 
charge density certainly influence the final baryon charge 
multiplicity and specific energy distributions. In the ab- 
sence of other dynamical effects during hadronization the 
fluctuations of the fluid dynamical densities will be inher- 
ited by the corresponding final measurable quantities. 

The present calculations show that in dynamical sys- 
tems, even with the least initial variation, strongly vary- 
ing higher statistical moments may develop. Thus the ef- 
fects of the FD expansion and of the final hadronization 
and freeze out should be separated. This can be done in 
theoretical hybrid models by evaluating separately both 
effects, and in experiments by mixed event methods (e.g. 
ref. |29() or more specific correlation measurements. 
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